In [412]:
def timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x, simulation=0,title=''):
plt.figure(figsize=(10, 6))
plt.subplot(311)
plt.plot(x, C[:,simulation], alpha = .5)
plt.ylabel('Coral Node Count')
plt.subplot(312)
plt.ylabel('CN Index')
plt.xlabel('Time')
plt.plot(x, CN[:,simulation], alpha = .5)
plt.subplot(313)
plt.ylabel('CNr Index')
plt.xlabel('Time')
plt.plot(x, CNr[:,simulation], alpha = .5)
def initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=1,
timesteps=[0,1,2,3,4,5,6,7,8,9,10,11,12]):
plt.figure(figsize=(15, 12))
#colors = ['pink', 'lightgreen','darkgreen']
#levels = [0, 1, 2]
#cmap, norm = clt.from_levels_and_colors(levels=levels, colors=colors, extend='max')
plt.subplot(331).set_title('-C- ' + str(C[timesteps[0],simulation])+
' -CN- ' + str(round(CN[timesteps[0],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[0],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[0],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(332).set_title('-C- ' + str(C[timesteps[1],simulation])+
' -CN- ' + str(round(CN[timesteps[1],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[1],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[1],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(333).set_title('-C- ' + str(C[timesteps[2],simulation])+
' -CN- ' + str(round(CN[timesteps[2],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[2],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[2],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(334).set_title('-C- ' + str(C[timesteps[3],simulation])+
' -CN- ' + str(round(CN[timesteps[3],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[0],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[3],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(335).set_title('-C- ' + str(C[timesteps[4],simulation])+
' -CN- ' + str(round(CN[timesteps[4],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[4],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[4],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(336).set_title('-C- ' + str(C[timesteps[5],simulation])+
' -CN- ' + str(round(CN[timesteps[5],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[5],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[5],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(337).set_title('-C- ' + str(C[timesteps[6],simulation])+
' -CN- ' + str(round(CN[timesteps[6],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[6],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[6],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(338).set_title('-C- ' + str(C[timesteps[7],simulation])+
' -CN- ' + str(round(CN[timesteps[7],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[7],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[7],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(339).set_title('-C- ' + str(C[timesteps[8],simulation])+
' -CN- ' + str(round(CN[timesteps[8],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[8],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[8],:,simulation],10),
cmap=cmap, norm=norm)
plt.show()
def initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=1,
timesteps=[0,1,2]):
plt.figure(figsize=(15, 4))
#colors = ['pink', 'lightgreen','darkgreen']
#levels = [0, 1, 2]
#cmap, norm = clt.from_levels_and_colors(levels=levels, colors=colors, extend='max')
plt.subplot(131).set_title('-C- ' + str(C[timesteps[0],simulation])+
' -CN- ' + str(round(CN[timesteps[0],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[0],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[0],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(132).set_title('-C- ' + str(C[timesteps[1],simulation])+
' -CN- ' + str(round(CN[timesteps[1],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[1],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[1],:,simulation],10),
cmap=cmap, norm=norm)
plt.subplot(133).set_title('-C- ' + str(C[timesteps[2],simulation])+
' -CN- ' + str(round(CN[timesteps[2],simulation],2))+
' -CNr- ' + str(round(CNr[timesteps[2],simulation],2)),
fontweight="bold")
plt.imshow(shaper(types[timesteps[2],:,simulation],10),
cmap=cmap, norm=norm)
plt.show()
Welcome to the Coral Model Testing and Exploration Notebook!¶
In this notebook, we take you through examples of running our model the various way to plot the model's outputs.
We begin with loading the necessary libraries¶
In [1]:
from coralModelTracking import Organism, Reef, Ocean
import tools as tl
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as clt
from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
Next, we set out model parameters¶
In [2]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .5
algaePercent = .25
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
Now we create and run multiple reefs¶
In [3]:
Pacific = Ocean()
for s in range(0,NumberOfSimulations):
Moorea = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([0,1,2],
p=[coralPercent, turfPercent, algaePercent])
node = Organism(type=U, location=[i,j], ID=count)
Moorea.append(node)
count = count + 1
Moorea.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
Moorea.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
Pacific.append(Moorea)
The Products we now can play with¶
Pacific.simulation[1].coralNodeCountPacific.simulation[1].coralNeighborCount
Functions to pull information from these objects¶
Can pull keys and values via:
x = list(dictionary.keys()), andy = list(dictionary.values())
In running the coralModel, we are interested in capturing the CN index: the average number of coral neighbors for each coral node.
Single Reef Example, looking at Simulation 1¶
In [4]:
example = Pacific.simulation[1].coralNodeCount
x = list(example.keys()) ## the time steps
coralCount = list(example.values()) ## the number of coral nodes
plt.figure(figsize=(15, 5))
plt.plot(x,coralCount)
Out[4]:
In [5]:
example = Pacific.simulation[1].coralNeighborCount
x = list(example.keys())
CN = np.array(list(example.values()))/np.array(coralCount) ## taking the average of the neighborhood count
plt.figure(figsize=(15, 5))
plt.plot(x,CN)
plt.plot(x,np.repeat(4, NumberOfTimesteps))
Out[5]:
Each node can have up to 8 neighbors. If the average is above 4 (the orange line above), this means on average, each coral node is surrounded by more coral than all other benthic types combined.
In [6]:
example = Pacific.simulation[1].coralNeighborCount
x = coralCount
y = CN
plt.figure(figsize=(15, 5))
plt.scatter(x,y)
Out[6]:
Looking at the Average across the various simulations¶
In [7]:
def dictToNumpy(dictionary):
output = list(dictionary.values())
output = np.array(output)
return(output)
multsimCoralCounts = np.array([dictToNumpy(Pacific.simulation[i].coralNodeCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(Pacific.simulation[i].coralNeighborCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
Coral Node Count: 7 simulations and average of 100 simulations¶
In [10]:
x = list(Pacific.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[10]:
CN index: 7 simulations and average of 100 simulation¶
In [11]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[11]:
Simulation with Maximized Coral Neighborhood Blob¶
In [12]:
## distance to center mask
center = (round(rows/2),round(columns/2))
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
plt.figure(figsize=(15, 5))
plt.imshow(distanceGrid)
Out[12]:
Checking Initial Reef¶
In [64]:
coralNodeLocations = (np.where(distanceGrid < np.median(distanceGrid)))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
### Sudo Making of Initial Reef
grid = np.zeros((rows,columns))
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
grid[i,j] = U
plt.figure(figsize=(15, 5))
plt.imshow(grid)
Out[64]:
Model set-up/run with initial coral blob condition, coral at 50%¶
In [15]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .5
algaePercent = .25
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Specific Coral Nodes
center = (round(rows/2),round(columns/2))
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
coralNodeLocations = (np.where(distanceGrid < np.median(distanceGrid)))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [16]:
## Running Model
Resilient = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
Resilient.append(optimalReef)
Coral Node Count: 7 simulations and average of 100 simulations¶
In [65]:
multsimCoralCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNodeCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNeighborCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
##plotting
x = list(Resilient.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[65]:
CN index: 7 simulations and average of 100 simulations¶
In [19]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[19]:
Model set-up/run with initial coral blob condition, coral at 33%¶
In [20]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]
coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [21]:
%%time
Resilient = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
Resilient.append(optimalReef)
Coral Node Count: 7 simulations and average of 100 simulations¶
In [23]:
multsimCoralCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNodeCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNeighborCount)
for i,val in enumerate(Pacific.simulation)]).transpose()
x = list(Resilient.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[23]:
CN index: 7 simulations and average of 100 simulations¶
In [24]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[24]:
G parameter = .5, Model set-up/run with initial coral blob condition , coral at 33%¶
In [80]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]
coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [81]:
%%time
LowGrazing = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
LowGrazing.append(optimalReef)
Coral Node Count: 100 simulations and average of 100 simulations¶
In [82]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[82]:
CN index: 100 simulations and average of 100 simulations¶
In [83]:
a = multsimNeighborCounts
b = multsimCoralCounts
CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)
plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[83]:
Same as above but with more time steps, testing with less simulations¶
In [84]:
## Simulations and Number of runs
NumberOfSimulations = 10
dt=.1
tf=1000
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]
coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [85]:
%%time
LowGrazing = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
LowGrazing.append(optimalReef)
In [86]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations¶
In [87]:
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[87]:
^^ For this one above it would be interesting to explore whether the spatial configuration for the two top lines allowed for coral success.¶
CN index: 10 simulations of 10,000 time steps, and the average of the 10 simulations¶
In [88]:
a = multsimNeighborCounts
b = multsimCoralCounts
CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)
plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[88]:
G parameter = .2, more time steps, testing with less simulations¶
In [89]:
## Simulations and Number of runs
NumberOfSimulations = 10
dt=.1
tf=1000
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.2
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]
coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [90]:
%%time
LowGrazing = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
LowGrazing.append(optimalReef)
In [91]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount)
for i,val in enumerate(LowGrazing.simulation)]).transpose()
Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations. Very grazing (g=.2)¶
In [92]:
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[92]:
Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations. Very low grazing (g=.2)¶
In [93]:
a = multsimNeighborCounts
b = multsimCoralCounts
CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)
plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[93]:
-- 30 SIMULATIONS, GRAZING .5, SPATIAL PLOTS ---¶
In [266]:
## Simulations and Number of runs
NumberOfSimulations = 30
dt=.1
tf=1000
## Grid shapes
rows, columns = 10, 10
## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33
## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9
threshold = 1.45 #extent of neighborhood recognition
NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent
## Storing Spatial Grid, may slow down run
#types = np.zeros((NumberOfTimesteps, NumberOfNodes, NumberOfSimulations))
## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center)
for i in range(0,rows)
for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]
coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n])
for n in range(0,len(coralNodeLocations[0]))]
In [267]:
%%time
#Switching = Ocean()
for s in range(0,NumberOfSimulations):
optimalReef = Reef()
count = 0
for i in range(0,rows):
for j in range(0,columns):
U = np.random.choice([1,2],
p=[.5, .5])
if (i,j) in coralNodeLocations:
U = 0
node = Organism(type=U, location=[i,j], ID=count)
optimalReef.append(node)
count = count + 1
optimalReef.generateGraph(threshold)
for n in range(0,NumberOfTimesteps):
for i,val in enumerate(optimalReef.nodes):
types[n,i,s] = optimalReef.nodes[i].type
optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt)
Switching.append(optimalReef)
In [274]:
# Save types, i.e. all the grid progression data
#for i in range(0,len(Switching.simulation)):
# np.savetxt("coralModelTrackingOutputs /modelOutput_switching"+ str(i)+".csv",
# np.reshape(types[:,:,i], (-1, rows)), delimiter=",")
In [441]:
#import pickle
#path = 'coralModelTrackingOutputs/'
#output_name = 'grazing5'
#outfile = open(path+output_name, "wb") #open pickle jar
#pickle.dump(Switching, outfile) #put contents into jar
#outfile.close() #close the jar
#PacificReturns = pickle.loads(open(path+output_name, "rb").read())
Coral Count over time, 30 Simulations¶
In [376]:
x = list(Switching.simulation[1].coralNodeCount.keys())
multsimCoralCounts = np.array([dictToNumpy(Switching.simulation[i].coralNodeCount)
for i,val in enumerate(Switching.simulation)]).transpose()
multsimNeighborCounts = np.array([dictToNumpy(Switching.simulation[i].coralNeighborCount)
for i,val in enumerate(Switching.simulation)]).transpose()
C = multsimCoralCounts
N = multsimNeighborCounts
plt.figure(figsize=(15, 5))
plt.plot(x,C, alpha=0.4)
plt.plot(x,C.mean(axis=1), 'black')
Out[376]:
CN index over time, 30 Simulations¶
In [351]:
CN = np.divide(N, C, out=np.zeros_like(N), where=C!=0)
CNmean = CN.mean(axis=1)
plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[351]:
CNr index over time, 30 Simulations¶
In [352]:
CNr = np.divide(CN, C, out=np.zeros_like(a), where=b!=0)
CNrmean = CNr.mean(axis=1)
plt.figure(figsize=(15, 5))
plt.plot(x, CNr, alpha=0.4)
plt.plot(x, CNrmean, 'black')
Out[352]:
Conclusions from exploring each of the simulations:¶
- Lives On (simulation:time at which relative steadyness is reached)
simulation 7:400, 14:1100, 25:900
Mid Crash (simulation:crashing time)
simulation 2:2300-3000, 22:3000, 28:3000
EarlyMid Crash (simulation:crashing time)
simulation 8:1800, 13:1700, 17:1800, 21, 23
Ambiguous early crashers (simulation:crashing time)
- simulation 0:1000, 1:300, 3:400, 4:300, 5:800, 6, 9, 10, 11, 12, 15, 16, 18, 19, 20, 24, 26, 27, 29
One task for tomorrow: run these grid configurations again and see if how their outputs vary (maybe this is where to include the variance)?¶
Plots where coral prevails - simulations 7, 14, 25¶
C, CN, CNr Time series example simulation 14¶
In [444]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=14,
title="Makes it for a while, and then crashes. grazing =.5")
Spatial evolutions of simulations 7, 14, and 25¶
In [464]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=7,
timesteps=[0,100,500])
In [465]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=14,
timesteps=[0,1000,2000])
In [466]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=25,
timesteps=[0,700,1000])
Plots from mid-crashes, simulations 2, 22, 28¶
C, CN, CNr Time series example simulation 8¶
In [417]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=2,
title="Makes it for a while, and then crashes. grazing =.5")
Spatial evolutions of simulations 2, 22, and 28¶
In [461]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=2,
timesteps=[0,2000,3000])
In [462]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=22,
timesteps=[0,2000,3000])
In [463]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=28,
timesteps=[0,2000,3000])
Plots from Early-Mid Crashes, simulations 8, 13, 17¶
C, CN, CNr time series example simulation 8¶
In [428]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=8,
title="Makes it for a while, and then crashes. grazing =.5")
Spatial evolutions of simulations 8, 13, 17¶
In [467]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=8,
timesteps=[0,1000,2000])
In [468]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=13,
timesteps=[0,1000,2000])
In [469]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=17,
timesteps=[0,1000,2000])
Plots from Early Crashes, simulations 6, 9, 10¶
C, CN, CNr time series example simulation 6¶
In [432]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=6,
title="Makes it for a while, and then crashes. grazing =.5")
Spatial evolutions of simulations 6, 9, 10¶
In [470]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=6,
timesteps=[0,100,300])
In [471]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=9,
timesteps=[0,100,300])
In [472]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=10,
timesteps=[0,100,300])
Looking at plots with same C but heading towards different states¶
Time Series simulation 14¶
In [487]:
timeSeriesPlotter(C=C[0:3000,:],CN=CN[0:3000,:],
CNr=CNr[0:3000,:],x=x[0:3000],simulation=14)
In [459]:
initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=14,
timesteps=[0,10,75,
90,110,125,
230,240,250])
Time Series simulation 2 (crashes)¶
In [488]:
timeSeriesPlotter(C=C[0:3000,:],CN=CN[0:3000,:],
CNr=CNr[0:3000,:],x=x[0:3000],simulation=2)
In [460]:
initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=2,
timesteps=[0,10,75,
90,110,125,
230,240,250])
Exploration of weird CNr vs CN effect¶
Coral Dominating Simulations¶
In [490]:
simulation = 7
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[490]:
In [491]:
simulation = 14
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[491]:
In [492]:
simulation = 25
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[492]:
Mid Crashes¶
In [493]:
simulation = 2
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[493]:
In [495]:
simulation = 22
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[495]:
In [496]:
simulation = 28
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[496]:
Coral Early-Mid Crashes¶
In [497]:
simulation = 8
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[497]:
In [498]:
simulation = 13
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[498]:
In [499]:
simulation = 17
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[499]:
Coral Early Crashes¶
In [501]:
simulation = 6
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[501]:
In [502]:
simulation = 9
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[502]:
In [503]:
simulation = 10
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[503]:
In [ ]: